1 intro

2 ‘ONs’

## network ----

d_subj_netw <- get_trial_data("network_ons") %>% get_ons
d_subj_netw_prw <- get_trial_data("network_ons_prw") %>% get_ons

# table(d_subj_netw$subj, is.na(d_subj_netw$value))


## plot:

d_subj_netw %>%
  
  ggplot(aes(roi, on)) +
  
  geom_hline(yintercept = 0) +
  stat_summary(aes(fill = task), fun = "mean", geom = "col", width = 0.5, position = position_dodge(width = 0.5)) +
  stat_summary(
    fun.data = "mean_cl_boot", geom = "errorbar", width = 0, size = 1.5, position = position_dodge(width = 0.5)
    ) +
  
  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none') +
  
  labs(y = "Mean ONs contrast (b)", x = "Network (Schaefer 7)", title = "vanilla")

d_subj_netw_prw %>%
  
  ggplot(aes(roi, on)) +
  
  geom_hline(yintercept = 0) +
  stat_summary(aes(fill = task), fun = "mean", geom = "col", width = 0.5, position = position_dodge(width = 0.5)) +
  stat_summary(
    fun.data = "mean_cl_boot", geom = "errorbar", width = 0, size = 1.5, position = position_dodge(width = 0.5)
    ) +
  
  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none') +
  
  labs(y = "Mean ONs contrast (b)", x = "Network (Schaefer 7)", title = "prewhitened")

d_subj_netw <- merge(d_subj_netw_prw, d_subj_netw, by = c("subj", "task", "roi"), suffix = c("_prw", "_van"))

d_subj_netw %>%
  
  melt(id.vars = c("subj", "task", "roi"), value.name = "on") %>%
  ggplot(aes(roi, on, fill = variable)) +
  
  geom_hline(yintercept = 2) +
  stat_summary(
    fun = function(x) t.test(x)$statistic, geom = "col", width = 0.5, position = position_dodge(width = 0.5)
    ) +
  
  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  # theme(legend.position = 'none') +
  
  labs(y = "Mean ONs contrast (t-stat)", x = "Network (Schaefer 7)")

## parcel ----

d_subj_parc_on <- get_trial_data("parcel_ons") %>% get_ons

## get parcel t-stats:

d_group_parc_on <- d_subj_parc_on[, .(tstat = t.test(on)$statistic, p = t.test(on)$p.value), by = c("task", "roi")]
d_group_parc_on[, p_fdr := p.adjust(p, "fdr"), by = task]  ## adjust p val
d_group_parc_on[, roi_ind := match(roi, key_schaefer$parcel)]  ## bind indices



d_subj_parc_on_prw <- get_trial_data("parcel_ons_prw") %>% get_ons

## get parcel t-stats:

d_group_parc_on_prw <- d_subj_parc_on_prw[, .(tstat = t.test(on)$statistic, p = t.test(on)$p.value), by = c("task", "roi")]
d_group_parc_on_prw[, p_fdr := p.adjust(p, "fdr"), by = task]  ## adjust p val
d_group_parc_on_prw[, roi_ind := match(roi, key_schaefer$parcel)]  ## bind indices

# any(is_equal(d_subj_parc_on_prw$on, d_subj_parc_on$on))

2.1 unthresholded

2.1.1 vanilla


for task_i in range(len(r.tasks)):
  # task_i = 0
  
  ## get task:
  
  is_task_i = r.d_group_parc_on.task == r.tasks[task_i]
  d_on = r.d_group_parc_on[is_task_i]
  
  ## get overlay: 

  overlay_on_lh = get_overlay(d_on.roi_ind, d_on.tstat, "left")
  overlay_on_rh = get_overlay(d_on.roi_ind, d_on.tstat, "right")
  overlay_on = np.column_stack(np.stack((overlay_on_lh, overlay_on_rh)))
  
  ## plot:
  
  fig = plot_surf_roi_montage(roi_map = overlay_on, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.

2.1.2 prewhitened


for task_i in range(len(r.tasks)):
  # task_i = 0
  
  ## get task:
  
  is_task_i = r.d_group_parc_on_prw.task == r.tasks[task_i]
  d_on = r.d_group_parc_on_prw[is_task_i]
  
  ## get overlay: 

  overlay_on_lh = get_overlay(d_on.roi_ind, d_on.tstat, "left")
  overlay_on_rh = get_overlay(d_on.roi_ind, d_on.tstat, "right")
  overlay_on = np.column_stack(np.stack((overlay_on_lh, overlay_on_rh)))
  
  ## plot:
  
  fig = plot_surf_roi_montage(roi_map = overlay_on, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
## <string>:30: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.

2.2 thresholded

2.2.1 vanilla


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_on.task == r.tasks[task_i]
  d_on = r.d_group_parc_on[is_task_i]

  ## get overlay:

  d_on.loc[d_on['p_fdr'] > 0.05, 'tstat'] = 0

  overlay_on_lh = get_overlay(d_on.roi_ind, d_on.tstat, "left")
  overlay_on_rh = get_overlay(d_on.roi_ind, d_on.tstat, "right")
  overlay_on = np.column_stack(np.stack((overlay_on_lh, overlay_on_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_on, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')

2.2.2 prewhitened


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_on_prw.task == r.tasks[task_i]
  d_on = r.d_group_parc_on_prw[is_task_i]

  ## get overlay:

  d_on.loc[d_on['p_fdr'] > 0.05, 'tstat'] = 0

  overlay_on_lh = get_overlay(d_on.roi_ind, d_on.tstat, "left")
  overlay_on_rh = get_overlay(d_on.roi_ind, d_on.tstat, "right")
  overlay_on = np.column_stack(np.stack((overlay_on_lh, overlay_on_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_on, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')

3 hi – lo

## network ----

d_subj_netw <- get_trial_data("network_hilo") %>% get_hilo
# table(d_subj_netw$subj, is.na(d_subj_netw$value))
d_subj_netw_prw <- get_trial_data("network_hilo_prw") %>% get_hilo

## plot:

d_subj_netw %>%

  ggplot(aes(roi, hi_v_lo)) +

  geom_hline(yintercept = 0) +
  stat_summary(aes(fill = task), fun = "mean", geom = "col", width = 0.5, position = position_dodge(width = 0.5)) +
  stat_summary(
    fun.data = "mean_cl_boot", geom = "errorbar", width = 0, size = 1.5, position = position_dodge(width = 0.5)
    ) +

  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none') +

  labs(y = "Mean hi vs lo contrast (b)", x = "Network (Schaefer 7)")

d_subj_netw_prw %>%
  
  ggplot(aes(roi, hi_v_lo)) +

  geom_hline(yintercept = 0) +
  stat_summary(aes(fill = task), fun = "mean", geom = "col", width = 0.5, position = position_dodge(width = 0.5)) +
  stat_summary(
    fun.data = "mean_cl_boot", geom = "errorbar", width = 0, size = 1.5, position = position_dodge(width = 0.5)
    ) +

  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none') +

  labs(y = "Mean hi vs lo contrast (b)", x = "Network (Schaefer 7)", title = "prewhitened")

d_subj_netw <- merge(d_subj_netw_prw, d_subj_netw, by = c("subj", "task", "roi"), suffix = c("_prw", "_van"))

d_subj_netw[, c("subj", "task", "roi", "hi_v_lo_prw", "hi_v_lo_van")] %>%
  
  melt(id.vars = c("subj", "task", "roi"), value.name = "hi_v_lo") %>%
  ggplot(aes(roi, hi_v_lo, fill = variable)) +
  
  geom_hline(yintercept = 2) +
  stat_summary(
    fun = function(x) t.test(x)$statistic, geom = "col", width = 0.5, position = position_dodge(width = 0.5)
    ) +
  
  facet_grid(vars(task), scales = "free_y") +
  scale_fill_brewer(type = "qual", palette = 2) +
  # theme(legend.position = 'none') +
  
  labs(y = "Mean hi vs lo contrast (t-stat)", x = "Network (Schaefer 7)")

# d_subj_netw[task %in% c("Axcpt", "Stroop"), c("subj", "task", "roi", "hi_v_lo_prw", "hi_v_lo_van")] %>%
#   # group_by(roi, subj) %>%
#   pivot_wider(id_cols = c("subj", "roi"), names_from = task, values_from = hi_v_lo_prw) %>%
#   group_by(roi) %>%
#   summarize(r = cor(Axcpt, Stroop))
  # summarize(r = cor(Cuedts, Stern))
## parcel ----

d_subj_parc_hilo <- get_trial_data("parcel_hilo") %>% get_hilo


## plot MD parcels:

d_subj_parc_hilo[roi %in% key_schaefer[schaefermd]$parcel] %>%

  ggplot(aes(task, hi_v_lo, fill = task)) +

  geom_hline(yintercept = 2) +
  stat_summary(
    fun = function(x) t.test(x)$statistic, geom = "col", width = 0.5, position = position_dodge(width = 0.5)
    ) +

  facet_wrap(vars(roi)) +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none', axis.text.x = element_blank()) +

  labs(y = "Mean hi vs lo contrast (t-stat)", title = "Schaeferized MMP-MD", x = "Task (alphabetical)")

## plot surfaces

## get parcel t-stats:

d_group_parc_hilo <- d_subj_parc_hilo[, .(tstat = t.test(hi_v_lo)$statistic, p = t.test(hi_v_lo)$p.value),
                                      by = c("task", "roi")]
d_group_parc_hilo[, p_fdr := p.adjust(p, "fdr"), by = task]  ## adjust p val
d_group_parc_hilo[, roi_ind := match(roi, key_schaefer$parcel)]  ## bind indices


## prewhitened ----

d_subj_parc_hilo_prw <- get_trial_data("parcel_hilo_prw") %>% get_hilo


## plot MD parcels:

d_subj_parc_hilo_prw[roi %in% key_schaefer[schaefermd]$parcel] %>%

  ggplot(aes(task, hi_v_lo, fill = task)) +

  geom_hline(yintercept = 2) +
  stat_summary(
    fun = function(x) t.test(x)$statistic, geom = "col", width = 0.5, position = position_dodge(width = 0.5)
    ) +

  facet_wrap(vars(roi)) +
  scale_fill_brewer(type = "qual", palette = 2) +
  theme(legend.position = 'none', axis.text.x = element_blank()) +

  labs(y = "Mean hi vs lo contrast (t-stat)", title = "Schaeferized MMP-MD | prewhitened", x = "Task (alphabetical)")

## plot surfaces

## get parcel t-stats:

d_group_parc_hilo_prw <- d_subj_parc_hilo_prw[, .(tstat = t.test(hi_v_lo)$statistic, p = t.test(hi_v_lo)$p.value),
                                      by = c("task", "roi")]
d_group_parc_hilo_prw[, p_fdr := p.adjust(p, "fdr"), by = task]  ## adjust p val
d_group_parc_hilo_prw[, roi_ind := match(roi, key_schaefer$parcel)]  ## bind indices

3.1 unthresholded

3.1.1 vanilla


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_hilo.task == r.tasks[task_i]
  d_hilo = r.d_group_parc_hilo[is_task_i]

  ## get overlay:

  overlay_hilo_lh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "left")
  overlay_hilo_rh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "right")
  overlay_hilo = np.column_stack(np.stack((overlay_hilo_lh, overlay_hilo_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_hilo, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')

3.1.2 prewhitened


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_hilo_prw.task == r.tasks[task_i]
  d_hilo = r.d_group_parc_hilo_prw[is_task_i]

  ## get overlay:

  overlay_hilo_lh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "left")
  overlay_hilo_rh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "right")
  overlay_hilo = np.column_stack(np.stack((overlay_hilo_lh, overlay_hilo_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_hilo, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')

3.2 thresholded

3.2.1 vanilla


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_hilo.task == r.tasks[task_i]
  d_hilo = r.d_group_parc_hilo[is_task_i]

  ## get overlay:

  d_hilo.loc[d_hilo['p_fdr'] > 0.05, 'tstat'] = 0

  overlay_hilo_lh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "left")
  overlay_hilo_rh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "right")
  overlay_hilo = np.column_stack(np.stack((overlay_hilo_lh, overlay_hilo_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_hilo, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')

3.2.2 prewhitened


for task_i in range(len(r.tasks)):
  # task_i = 0

  ## get task:

  is_task_i = r.d_group_parc_hilo_prw.task == r.tasks[task_i]
  d_hilo = r.d_group_parc_hilo_prw[is_task_i]

  ## get overlay:

  d_hilo.loc[d_hilo['p_fdr'] > 0.05, 'tstat'] = 0

  overlay_hilo_lh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "left")
  overlay_hilo_rh = get_overlay(d_hilo.roi_ind, d_hilo.tstat, "right")
  overlay_hilo = np.column_stack(np.stack((overlay_hilo_lh, overlay_hilo_rh)))

  ## plot:

  plot_surf_roi_montage(roi_map = overlay_hilo, title = r.tasks[task_i])
  plotting.show()
  fig.clear()
  plt.close('all')